Vector velocity estimation using directional beam forming and cross correlation

ABSTRACT

The two-dimensional velocity vector using a pulsed ultrasound field can be determined with the invention. The method uses a focused ultrasound field along the velocity direction for probing the moving medium under investigation. Several pulses are emitted and the focused received fields along the velocity direction are cross-correlated. The time shift between received signals is found from the peak in the cross-correlation function and the velocity is thereby determined.

FIELD OF THE INVENTION

The invention relates to an apparatus and a method for determining thevelocity vector of a remotely sensed object using either sound, inparticular ultrasound, or electromagnetic radiation. The movement of theobject is determined by emitting and receiving a pulsed field focused inthe direction of the velocity vector. By using a number of pulseemissions, the inter pulse movement can be estimated and the velocityfound from the estimated movement and the time between pulses. Theinvention is based on the principle of using a directionally focusedfield for making the received signal influenced by motion along thedirection of the movement.

BACKGROUND OF THE INVENTION

Medical ultrasound is extensively used for studying flow dynamics in thehuman body by using color flow mapping. The technique displays a colorimage of the flow superimposed on the normal anatomic B-mode image.Traditionally, the velocity component along the ultrasound beamdirection is measured, and a flow transverse to the beam is notdisplayed, since it is not measured or estimated. An example of this isshown in FIG. 1, where the flows in the carotid artery and the jugularvein are displayed. The image is acquired with a convex arraytransducer, and the angles between flow direction and the ultrasoundbeam change over the image. Notice the change of estimated flowdirection around the dashed line in both vessels due to the change ofangle between the flow and the ultrasound beam. This is one of the mainlimitations of current ultrasound flow systems, since most vessels areparallel to the skin surface, and therefore it is a problem to get asufficiently small angle between the flow and the beam. Also the flow isoften not parallel to the vessel surface, and it is therefore difficult,if not impossible, to estimate the correct angle and compensate for it[1]. In European patent application EP 616 231 [2] the velocity is foundthrough a cross-sectional area using a 2D matrix transducer that canfocus on the individual areas in the cross-section. The volume flowthrough the cross-section is then found, but still only the velocity indirection of the ultrasound beam is estimated. Several authors haveattempted to remedy this artifact. Fox [3] suggested using two beams tofind the transverse component. The system works well for largetransducers and investigations close to the transducer, but the varianceof the transverse component increases for situations with large depthsand smaller transducers as used in cardiac scanning through the ribs.Trahey and co-workers [4] have suggested using speckle tracking in whicha small search region in one image is correlated or compared to asubsequent image. This approach has problems in terms of frame rate,since images are compared, and the resolution of the velocity estimatescan be low. Newhouse et al. [5] developed a method in which the totalbandwidth of the received signal is affected by the transverse velocity.It is, however, often difficult to find this bandwidth due to theinherent noise in the signal.

Of special interest is the working by Bonnefous [6], which uses a numberof parallel beams to find the transverse velocity. The approach does,however, not work for a velocity that is not orthogonal to theultrasound beam direction.

In this invention a new technique using a focused signal in thedirection of the velocity is used. The velocity is then found byacquiring two or more of these focused signals and cross correlatingthem to find the displacement between pulse emission, whereby thevelocity can be determined.

PRIOR ART APPROACH

This section summarizes the article written in 1988 by Bonnefous [6],where a method for estimating the transverse velocity was suggested.

Transverse velocity estimation must perform a signal processing, wherethe effect of axial motion is negligible compared to the transverse one.The idea presented by Bonnefous requires a broad beam in emission (aplane ultrasound wave front), and a number of parallel identical beams,separated by a pitch w in the transverse direction, are generated inreception, see FIG. 2.

For a given depth z, the signal received S_(n)(x,t) from the beamcentered in x=0 is

 S _(n)(0,t)=∫p(x,t)D _(n)(x)dx  (1)

where p(x,t) is the response of a scatterer located at x in theultrasound beam centered at x=0, and D^(n)(x) is the scattererdistribution at the instant nT_(prf), where T_(prf) is the pulserepetition period. In the same way, the signal received from the beamcentered at x=w is

S _(n)(w,t)=∫p(x−w,t)D _(n)(x)dx  (2)

If only a transverse uniform motion of the scatterers in considered, thedisplacement between two consecutive pulses nT_(prf) and (n+1)T_(prf)will give the relation for the distribution of the scatterers

D _(n+1)(x)=D _(n)(x−v _(x) T _(prf))  (3)

where v_(x) is the velocity of the scatterers. Combining equations (2)and (3) gives

S _(n+1)(0,t)=∫p(x,t)D _(n+1)(x)dx=∫p(x,t)D _(n)(x−v _(x) T_(prf))dx=∫p(x+v _(x) T _(prf) ,t)D _(n)(x)dxS _(n+1)(0,t)=S _(n)(−v_(x) T _(prf) ,t)  (4)

For a beam centered in x=w the relation between the consecutive signalsis

S _(n+1)(w,t)=S _(n)(w−v _(x) T _(prf) ,t)  (5)

The signal received in x=w at the pulse time (n+1)T_(prf) is, thus, thesame as the signal received in a beam centered in x=w−v_(x)T_(prf) atthe pulse time nT_(prf).

The correlation between the two signals p(x,t) and p(x−w,t) is anaveraging of the received signals over the random scattererdistributions. The cross-correlation of the received signals from twoadjacent signals is $\begin{matrix}\begin{matrix}{{C_{1}(w)} = {\sum\limits_{n}{\int_{t_{0}}^{t_{0} + {\Delta \quad t}}{{S_{n}\left( {0,t} \right)}{S_{n + 1}\left( {w,t} \right)}\quad {t}}}}} \\{= {\sum\limits_{n}{\int_{t_{0}}^{t_{0} + {\Delta \quad t}}{{S_{n}\left( {0,t} \right)}{S_{n}\left( {{w - {v_{x}T_{prf}}},t} \right)}\quad {t}}}}} \\{{C_{1}(w)} = {C_{0}\left( {{w - {v_{x}T_{prf}}},t} \right)}}\end{matrix} & (6)\end{matrix}$

where C₀(w) = ∑∫_(t₀)^(t₀ + Δ  t)S_(n)(0, t)S_(n)(w, t)  t,

is the autocorrelation function averaged over a number of receivedlines, where the line number is denoted by n. The interval (t₀,t₀+Δt) isthe range gate selected for the received signals. Equation (6) show thatthe shift of C₁ compared with C₀ is the transverse displacement betweenthe instants nT_(prf) and (n+1)T_(prf). Therefore, the maximum of C₁(w)is C₁(v_(x)T_(prf))=C₀(0).

For the general case, in which both axial and transverse motion takesplace, the equation that relates the received signals from twosuccessive pulses will be $\begin{matrix}{{S_{n + 1}\left( {w,t} \right)} = {S_{n}\left( {{w - {v_{x}T_{prf}}},{t - \quad \frac{2v_{z}T_{prf}}{c}}} \right)}} & (7)\end{matrix}$

Here t_(s)=2v_(z)T_(prf)/c, is the time shift for the axial motion.

The cross-correlation and autocorrelation functions are generalized totwo-dimensional functions and their expressions are $\begin{matrix}\begin{matrix}\begin{matrix}{{C_{1}\left( {w,u} \right)} = \quad {\sum\limits_{n}{\int_{t_{0}}^{t_{0 + {\Delta \quad t}}}{{S_{n}\left( {0,t} \right)}{S_{n + 1}\left( {w,{t + u}} \right)}\quad {t}}}}} & {\quad (8)} \\{= \quad {\sum\limits_{n}{\int_{t_{0}}^{t_{0} + {\Delta \quad t}}{{S_{n}\left( {0,t} \right)}{S_{n}\left( {{w - {v_{x}T_{prf}}},{t - \quad \frac{2v_{z}T_{prf}}{c} + u}} \right)}\quad {t}}}}} & {\quad (9)} \\{{C_{0}\left( {w,u} \right)} = \quad {\sum\limits_{n}{\int_{t_{0}}^{t_{0} + {\Delta \quad t}}{{S_{n}\left( {0,t} \right)}{S_{n}\left( {w,{t + u}} \right)}\quad {t}}}}} & {\quad (10)}\end{matrix} & \quad\end{matrix} & \quad\end{matrix}$

The relation between the cross-correlation and the autocorrelation is$\begin{matrix}{{C_{1}\left( {w,u} \right)} = {C_{0}\left( {{w - {v_{x}T_{prf}}},{u - \quad \frac{2v_{z}T_{prf}}{c}}} \right)}} & (11)\end{matrix}$

The two-dimensional determination of the velocity vector is performed byfirst finding the axial velocity and then the transverse component.

1. Axial velocity measurement: C₁(0,u) is calculated and then, the timeposition of the correlation peak is determined as $\begin{matrix}{{\underset{u}{Max}\left\lbrack {C_{1}\left( {0,u} \right)} \right\rbrack} = {C_{1}\left( {0,\quad \frac{2v_{z}T_{prf}}{c}} \right)}} & (12)\end{matrix}$

2. Transverse velocity measurement: Fixing the value of the timecoordinate to be u=2v_(z)T_(prf)/c in C₁(w,2v_(x)T_(prf)/c), the peakfor the transverse correlation is determined as $\begin{matrix}{\underset{w}{Max}\left\lbrack {{C_{1}\left( {w,\quad \frac{2v_{z}T_{prf}}{c}} \right\rbrack} = {C_{1}\left( {{v_{x}T_{prf}},\quad \frac{2v_{z}T_{prf}}{c}} \right)}} \right.} & (13)\end{matrix}$

This process can be generalized to measure the three-dimension velocityvector. New parallel beams along the y-axis should be used, and thecross-correlation function will be a three-dimensional function of thetype C₁(w,h,u). The problem with this approach is, however, that aspatially invariant velocity field is assumed, which is not the case inthe human body.

SUMMARY OF THE INVENTION

This section describes the approach of the invention for finding thetransverse component of the blood velocity. The invention presupposesthat the direction of movement is known. In medical ultrasound atwo-dimensional image including a blood vessel of interest is produced.The operator will then manually indicate the position and direction ofthe blood vessel to the system, or this information can be obtainedautomatically.

In a first embodiment of the invention, wave energy is emitted in apredetermined direction towards a moving object or a collection ofmoving objects, which will interact with the wave energy, whereby thewave energy will be scattered or reflected. Moving objects thusinteracting with emitted wave energy are often referred to as“scatterers”. Scattered or reflected wave energy will be received andprocessed to yield the desired estimate or measurement of velocity ofthe moving object or a collection of moving objects. The velocity of anyscattering or reflecting object, whether emitting it self or not, can bemeasured.

The main idea of the first embodiment of the invention is the creationof a plurality of focal points, which together form a focused line orbeam in the direction of the velocity vector that will track the motionof the scatterers. The generation of this focused beam requires a broadbeam in emission and multiple foci along the line in reception. This canbe achieved with a multi-element transducer. The focusing line issituated over the region of interest, where the motion of the scatterersoccurs, see FIG. 3. The method of the invention presupposes thatscatterers, whose motion is tracked, have the same velocity along thewhole length of the line. When measuring a laminar parabolic flow orother spatially variant velocity fields, the focus line therefore has tobe oriented in the direction of the flow lines.

The lateral beam is calculated during reception by delaying and addingresponses from the individual elements of the array. The delays arefound for each focus point by calculating the time it takes for theultrasound wave to travel from each transducer elements to the focuspoint and back to the transducer. Lateral beams are in this wayconstructed for each of the emitted pulses.

The velocity is estimated from the cross-correlation of two consecutivelateral signals. Hereby the displacement of the signals corresponds toan estimate of the distance traveled by the scatterers in the lateralbeam direction. Since the lateral beams are situated along the flowstream-lines, the complete velocity magnitude can be directly estimated.

In a second embodiment the velocity of a moving object or a collectionof moving objects emitting wave energy while moving, is measured. Noemission of wave energy is required, but rather the wave energy emittedby the moving object is received and processed to yield the desiredestimate or measurement of velocity of the moving object or a collectionof moving objects. The main idea of the second embodiment of theinvention is the creation of a plurality of focal points, which togetherform a focused line or beam in the direction of the velocity vector thatwill track the motion of the emitting object or objects.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow image of the cartiod artery and the jugular veinscanned with a convex array transducer;

FIG. 1A is a schematic block diagram of an apparatus in accordance withthe invention;

FIG. 2 illustrates parallel ultrasound beams in reception;

FIG. 3 illustrates focusing lines for obtaining lateral beams in thecase of a parabolic velocity profile;

FIG. 4 shows a coordinate system for the lateral beam measurementprinciple;

FIG. 5 shows delay graphs for the first, middle and last focus point ofa transverse line og length 10 mm perpendicular to the z-axis. Thedashed curves represent the delay lines compensated for the axialvelocity of the scatterers when θ=30°, v=0.5 m/s and T_(prf)=1/3500 s;

In FIG. 6 the top graph shows two consecutive lateral signals, where thesecond signal is compensated for the axial movement. The displacement,thus, corresponds to the transverse motion. In the bottom graph it hasbeen evaluated by the cross-correlation method. Simulation was performedfor θ=60°, v=0.5 m/s and dx=0.01 mm;

FIG. 7 shows the definition of depth in the vessel;

FIG. 8 shows the evolution of the lateral signals from two ultrasoundpulses: pulse (i) is the continuous line, and pulse (I+1) is the dashedline. Five different points of the cross-section of the vessel areshown. The distance is measured from the center axis of the vessel. Thedata correspond to an inclination of the vessel of θ=60°, see also FIG.12;

In FIG. 9 the continuous line arrows show the real velocity, and thedashed line arrows show the estimated velocity: v_(est)=({circumflexover (v)}_(t),{circumflex over (v)}_(z))for plug flow. The angle θ ismeasured clockwise from the vertical axis;

FIG. 10 shows estimated velocity profiles for θ=30° using 100, 20, 10and 4 lines. The thick lines indicate the true velocity profile, thecontinuous line is the mean of the estimates, and the dotted lines arethe standard deviations;

FIG. 11 shows estimated velocity profiles for θ=45° using 100, 20, 10and 4 lines. The thick lines indicate the true velocity profile, thecontinuous line is the mean of the estimates, and the dotted lines arethe standard deviations;

FIG. 12 shows estimated velocity profiles for θ=60° using 100, 20, 10and 4 lines. The thick lines indicate the true velocity profile, thecontinuous line is the mean of the estimates, and the dotted lines arethe standard deviations; and

FIG. 13 shows estimated velocity profiles for θ=90° using 100, 20, 10and 4 lines. The thick lines indicate the true velocity profile, thecontinuous line is the mean of the estimates, and the dotted lines arethe standard deviations.

DETAILED DESCRIPTION OF THE INVENTION

In FIG. 1 is shown an ultrasound flow image of the cartiod artery andthe jugular vein of a human being scanned with a convex arraytransducer. Flow towards the transducer and flow away from thetransducer will normally be shown in different colors such as red andblue, but here the colors are reproduced as different shades of gray.The dashed line is perpendicular to both blood vessels, and on bothsides of the dashed line a black region indicates that no flow isregistered in the direction of the dashed line. On both sides of theblack no-flow region different shades of gray indicate the direction offlow as marked with arrows in FIG. 1.

In FIG. 1A an apparatus 10 in accordance with the invention is shownschematically. A generator and pulser block 11 generates pulses ofultrasound frequency excitation signals. A transducer element selectorblock 12 receives the pulses of ultrasound frequency excitation signalsfrom the generator and pulser block 11 and feeds the signals to selectedones of the transducer elements of an emitting transducer 13. Theemitting transducer 13 is thereby caused to emit a beam 14 of ultrasoundsignals in a predetermined direction toward an object 15 underinvestigation. The object 15 under investigation interacts with the beam14 of ultrasound signals and scatters the beam. A portion of thescattered ultrasound field 15 is received in a receiving transducer 16.In the receiving transducer 16 transducer elements output electricsignals in dependence on the received ultrasound energy, and a beamformer 17 receives the electric signals and performs a beam forming. Thebeam formed signals from the beam former 17 are received in a processor18 that processes the beam formed signals into images 19.

The object 15 under investigation can be the body of an animal or humanbeing with blood vessel 20 of radius R in the body. In the blood vessel20 the blood flow has a velocity profile 21, e.g. as indicated in FIG.3. A plurality of focal points is combined to define a focusing line 22or focusing beam in then predetermined direction of flow as in FIG. 4.The focusing line 22 can be placed anywhere in the blood vessel 20, anda plurality of lateral beams 22 can be formed as in FIG. 3.

Measurement Principle of the Invention

This section derives a mathematical model of the tilted beam approach.The measurement situation is depicted in FIG. 4. The figure shows amulti-element transducer with N elements. The vector from the origin ofreference to the physical center of element i is {right arrow over(r)}_(i). The vector indicating the position of the focus points in theline is {right arrow over (r)}_(j) and is given by

{right arrow over (r)}_(i)=(x _(i) ,y _(i) ,z _(i)), i=1. . . N {rightarrow over (r)} _(j)=(x _(j) ,y _(j) ,x _(j)), j=1. . . M {right arrowover (r)} _(ij) ={right arrow over (r)} _(j) −{right arrow over (r)}_(i)=(x _(i) −x _(i) ,y _(j) −y _(i) ,z _(j) −z _(i)),  (14)

where {right arrow over (r)}_(ij) is the vector from element i to thefocus point j, and M is the number of focus points along the line. Thetime it takes for the ultrasound pulse to get from the element i to thefocus point j and back to the element i is

t _(ij)2/c|{right arrow over (r)} _(ij)|=2/c{square root over ((x _(j)−x _(i))²+(y _(j) −y _(i))²+(z _(j) −z _(i))²)}  (15)

The received backscattered signal in element i that corresponds to theultrasound pulse emitted in nT_(prf), can be written as r_(n)({rightarrow over (r)}_(i),t), where t indicates the time since pulse emission.

Assuming linearity the value of the field in a focus point j, iscalculated by adding the responses found in each element i for theappropriated delay time t_(ij) $\begin{matrix}{{s_{n}\left( {\overset{\rightarrow}{r}}_{j} \right)} = {\sum\limits_{i = 1}^{N}{r_{n}\left( {{\overset{\rightarrow}{r}}_{i},t_{ij}} \right)}}} & (16)\end{matrix}$

For the next ultrasound pulse (n+1)T_(prf), the value of the field inthe same point will be $\begin{matrix}\begin{matrix}{{s_{n + 1}\left( {\overset{\rightarrow}{r}}_{j} \right)} = {\sum\limits_{i = 1}^{N}{r_{n + 1}\left( {{\overset{\rightarrow}{r}}_{i},t_{ij}} \right)}}} \\{{= {\sum\limits_{i = 1}^{N}{r_{n}\left( {{\overset{\rightarrow}{r}}_{i},{t_{ij} - {2\quad \frac{v_{i,{axis}}}{T_{prf}}}}} \right)}}},}\end{matrix} & (17)\end{matrix}$

where v_(i,axis) represents the component of the velocity of a scattererplaced in {right arrow over (r)}_(j), projected on the {right arrow over(r)}_(ij) direction. The expression above uses the relation between twoconsecutive backscattered signals. The relationt_(ij)−2v_(i,axis)/T_(prf) will correspond to a delay value t_(ik),where k is another point on the focusing line. If{right arrow over(d)}_(l) is a unitary vector in the direction of the lateral beam, andl_(s) is the distance between the points j and k, then {right arrow over(r)}_(j={right arrow over (r)}) _(k)+l_(s){right arrow over (d)}_(l) andthe following relation is found $\begin{matrix}\begin{matrix}{{s_{n + 1}\left( {\overset{\rightarrow}{r}}_{j} \right)} = {\sum\limits_{i = 1}^{N}{r_{n}\left( {{\overset{\rightarrow}{r}}_{i},{t_{ij} - {2\quad \frac{v_{i,{axis}}}{T_{prf}}}}} \right)}}} \\{= {\sum\limits_{i = 1}^{N}{r_{n}\left( {{\overset{\rightarrow}{r}}_{i},t_{ik}} \right)}}} \\{{s_{n + 1}\left( {\overset{\rightarrow}{r}}_{j} \right)} = {{s_{n}\left( {\overset{\rightarrow}{r}}_{k} \right)} = {s_{n}\left( {{\overset{\rightarrow}{r}}_{j} - {l_{s}{\overset{\rightarrow}{d}}_{l}}} \right)}}}\end{matrix} & (18)\end{matrix}$

The beam processing is performed on digital signals and the sampledversion of the beam formed signal is denoted R_(n)[l]. The signalcalculated for the pulse at time (n+1)T_(prf) is a shifted version ofthe lateral signal for the pulse at time nT_(prf). The shift in positioncan be calculated from the cross-correlation of two consecutive signalsS₁[l] and S₂[l]: $\begin{matrix}\begin{matrix}{R_{12} = {\frac{1}{M}{\sum\limits_{l = 0}^{M - 1}{{S_{1}\lbrack l\rbrack}{S_{2}\left\lbrack {l + n} \right\rbrack}}}}} \\{= {\frac{1}{M}{\sum\limits_{l = 0}^{M - 1}{{S_{1}\lbrack l\rbrack}{S_{1}\left\lbrack {l + n - l_{s}} \right\rbrack}}}}} \\{= {R_{11}\left( {k - l_{s}} \right)}}\end{matrix} & (19)\end{matrix}$

The shift in position is then determined by the index of the peak of thecross-correlation function $\begin{matrix}{{\underset{k}{Max}\left\lbrack {R_{12}\quad (k)} \right\rbrack} = {\left. {\underset{k}{Max}\left\lbrack {R_{11}\quad \left( {k - l_{s}} \right)} \right\rbrack}\Rightarrow k \right. = {\hat{l}}_{s}}} & (20)\end{matrix}$

The cross-correlation can be improved by averaging over severalestimates of R₁₂, since the velocity of the scatterers can be consideredconstant for several pulses. The estimated mean velocity between twoconsecutive pulses is: $\begin{matrix}{\hat{v} = \frac{{\hat{l}}_{s}{dl}}{T_{prf}}} & (21)\end{matrix}$

where dl is the sampling interval along the lateral beam direction.

Functionality of the Invention

The invention was simulated using the simulating program Field II [7].The lateral beam is directly calculated with a focusing strategy thatcan be achieved with a regular array transducer. The method was firsttested for a blunt profile of the scatterers and then for a parabolicprofile.

Estimation of the Velocity for a Blunt Profile

In the beginning, the lateral beam was intended as an attempt toestimate the transverse velocity component ({circumflex over (v)}_(x)).Since {circumflex over (v)}_(z) could be obtained from the time-shiftestimation method, the velocity in the x-z plane will be estimated as{right arrow over (v)}_(est)=({circumflex over (v)}_(x),{circumflex over(v)}_(z)).

The estimation of v_(x) can be obtained in the case of a blunt profileusing a focusing line perpendicular to the z-axis. Two lateral signalsare generated for each ultrasound pulse. The second one has beencompensated for the axial movement that took place in the pulse intervaltime. The cross-correlation of the first signal obtained from pulse iwith the compensated signal from pulse i+1 reveals the transverse motionof the scatterers between the two pulses.

The calculation of the lateral beam and the estimation of the velocityare explained below. The last part shows the results obtained fordifferent vessel inclinations and for two different series ofparameters.

Simulations for a Blunt Profile

The simulation is performed in the following way:

1. The transducer used for the simulation is a linear array transducerwith 64 elements of dimensions: 0.25 mm width, 5 mm height and a kerf orseparation between elements of 0.05 mm. It gives the possibility ofelectronic focusing and apodization. A Hanning apodization is used forboth emission and reception. The emit focus is set to be at an infinitedistance from the transducer, so that a plane wave front is found at themeasurement area. Short ultrasound pulses are simulated to be emittedevery T_(prf) seconds.

2. The generation of simulated data must include an artificial phantomconsisting of point scatterers. This phantom simulates a blood vesselclose to the skin surface at a depth of 30 mm. The angle between theblood vessel and the axial ultrasound beam is called 0, see FIG. 4. Thephantom will be propagated a distance Δr=vT_(prf) in the intervalbetween two consecutive ultrasound pulse emissions to simulate thedisplacement of the scatterers. In the first simulations a blunt profileis considered, so that all scatterers move at the same velocity that wasset to be 0.5 m/s.

3. For every ultrasound pulse, an individual simulation for each elementof the transducer is done by the Field II program and stored. The fieldprogram is called 64 times for each of the transducer elements and foreach pulse emission.

4. The focused beam is calculated as a transverse beam situated at adepth of 30 mm, in the middle of the region of interest. It uses thedata obtained in the previous step. The line, transverse to theultrasound beam has the same length as the transducer, i.e. 20 mm. Theseparation between two focus points was chosen to be dx=0.02 mm, whichgives a total of 1000 samples for the transverse beam. The value of thefield in each focus point is found by calculating the time that it takesthe ultrasound wave to travel from each of the transducer elements tothe field point and back to the transducer. This is

delay=2/c(x _(f) −x _(c))²+(y _(f) −yc)²+(z _(f) −z _(c))²  (22)

where (x_(c),c_(y),z_(c)) are the positions of the center of thephysical elements of the aperture and (x_(f),y_(f),z_(f)) the positionof each point in the transverse line. Three delay lines for the first,middle and last points of the transverse line have been plotted in FIG.5 (continuous line).

The field value in each focus or transverse sampling point is the sum ofthe fields values obtained for each channel at the times thatcorresponds to the calculated delay line. The signal constructed thisway is stored for each ultrasound emission.

At the same time and for each emission, a second signal is synthesized.It corresponds to the lateral beam where the axial movement of thescatterers has been compensated for. In this case the delays are$\begin{matrix}{{delay}^{\prime} = {{\frac{2}{c}\sqrt{\left( {x_{f} - x_{c}} \right)^{2} + \left( {y_{f} - y_{c}} \right)^{2} + \left( {z_{f} - z_{c}} \right)^{2}}} - \quad {\frac{2v_{z}}{c}T_{prf}}}} & (23)\end{matrix}$

These new delay lines are plotted as dashed lines in FIG. 5. The axialvelocity (v_(z)) needs to be estimated before the calculation of thecompensated delay lines.

Estimation of the Velocity

The signals obtained from the simulation are used to estimate thevelocity. The velocities are estimated by cross-correlating twoconsecutive signals in which the second signal is compensated for theaxial motion that takes place during the pulse time interval. Theposition of the cross-correlation peak then indicates the shift of thesignals.

The transverse velocity is found from $\begin{matrix}{{{\hat{v}}_{x} = {\frac{dl}{T_{prf}} \cdot n_{int}}},} & (24)\end{matrix}$

where n_(int) is the interpolated index found by fitting a parabolaaround the maximum in the cross-correlation function [8].

Simulations and Results

FIG. 9 shows the result of the simulations. They were performed for ablunt profile and for a scatterer velocity equal to v=0.5 m/s. Thescatterer patterns were generated for seven different angles between thevessel and the axial beam. The received signals were found for 20different pulse emissions. The cross-correlation estimator used 4consecutive lines. This gave a total of 16 estimates for each velocityvalue. The dashed arrows point to (mean({circumflex over(v)}_(x)),mean({circumflex over (v)}_(z))) for each angle, while thecontinuous line arrows represent the true values. The standard deviationis represented for each estimated point as the axis of the followingellipse: $\begin{matrix}{{\left( \frac{{mean}\left( {\hat{v}}_{x} \right)}{{std}\left( {\hat{v}}_{x} \right)} \right)^{2} + \left( \frac{{mean}\left( {\hat{v}}_{z} \right)}{{std}\left( {\hat{v}}_{z} \right)} \right)^{2}} = 1} & (25)\end{matrix}$

The simulation parameters are shown in Table 1 and the results areplotted in FIG. 9 and listed in Table 2.

Estimation of the Velocity for a Parabolic Profile

The received signals are now focused along the flow direction. Thelateral beams are tilted in order to follow the motion of thescatterers, see FIG. 3. The inclination of the vessel has to be known orestimated with a B-mode image for the calculation of the flow beams.

In the simulation a box of 20×10×20 mm³ is created. It contains 10point-scatterers per cubic wavelength. The position and amplitude of thescatterers are randomly generated. The box of scatterers is rotated anangle φ=π/2−θ to simulate the vessel inclination. The ‘vessel’ has aradius of 5 mm and is situated at a depth of 30 mm from the transducer.The scatterers that lie within the vessel wall are propagated betweentwo consecutive pulses. The velocity profile is parabolic and thevelocity in the center of the vessel is 0.5 m/s.

For each ultrasound pulse i, the data, with the positions and amplitudesof the scatterers, is loaded. Then, a scan for each of the channels ofthe array transducer is done. Finally, the delay lines for each focusingline are calculated and the flow beams are generated by adding theresponses that correspond to the delay lines.

The transducer used was a linear array with 64 elements as described inSection 5.1. The parameters used for the simulation are listed in Table3.

The focused lines had a length of 10 mm. They are rotated to have thesame direction as the flow lines (see FIG. 7). They have to cover anaxial distance of at least $\begin{matrix}{{Z(\Theta)} = \frac{2R}{\sin \quad \Theta}} & (26)\end{matrix}$

This length assures that there will be flow-beams in all thecross-section of the vessel. For small values of θ, the axial range-gatehas to be increased. Z(θ) is divided in small segments, and for each ofthem a focusing line is calculated. The axial and lateral intervalvalues used in the simulation are given in Table 4.

The velocity estimates are found from the simulated data. The velocitymap is calculated by the cross-correlation of two received signals.These signals have been plotted in FIG. 8 for five different depths. Thecontinuous lines represents the first signals and the dashed lines thesecond. The motion of the scatterers can be appreciated by thedisplacement of the signals. In the center of the vessel (r=0) thelargest shift in position is observed. This shift decreases with thedistance to the axis, and a parabolic profile can be recognized.

Once the cross-correlation function is calculated and averaged over anumber of lines (4, 10, 20 or 100 for this simulation), the index of themaximum in the cross-correlation is obtained and interpolated, and thevelocity is found from (24).

The minimum velocity that can be obtained is limited by the lateralsampling interval (dl), and the maximum value is determined by the lagin the cross-correlation function: $\begin{matrix}{v_{\min} = \frac{dl}{T_{prf}}} & (27) \\{v_{\max} = {\frac{dl}{T_{prf}} \cdot {lag}}} & (28)\end{matrix}$

For the values introduced in the simulation: dl=0.01 mm,T_(prf=){fraction (1/3500)} s and a sampling transverse line of 10 mm,this gives 10/0.01=1000 sampling points and the limits for the velocityestimation are

v _(min)=0.035 m/s  (29)

v _(max)=0.035·1000=35 m/s  (30)

The simulation was repeated for four inclinations of the vessel, θ=30,45, 60, and 90°. For each angle, the estimates were found when averagingover 4, 10, 20 or 100 lines in the cross-correlation, which gave a totalof 96, 90, 80 or 1 velocity estimates respectively for every lateralbeam. The results are plotted in figures: 10, 11, 12, and 13. The thickline is the true velocity profile, the continuous line is the mean ofthe estimates, and the dotted lines are the standard deviations of theestimates.

The error in the velocity estimation is represented by the mean and thestandard deviation of the vertical distance between the true and theestimated velocity profile (Table 5). The values in the table are whenusing 4 lines for calculating the cross-correlation.

TABLE 1 Parameters used for the blunt profile simulation. The emit focaldistance is 150 mm. Pulse Number Center Wave- Number repeti- of Trans-fre- length of Sampling tion scat- versal quency λ = cycles frequencyperiod terers interval f₀[Hz] c/f₀ [m] M f_(s) [Hz] T_(prf) [s] N dx [m]3 10⁶ 0.51 10⁻³ 2 100 10⁶ 1/3500 5000 0.01 10⁻³

TABLE 2 Θ[deg] 0 30 60 90 135 −135 −45 True v_(x) [m/s] 0 0.25 0.43 0.50.35 −0.35 −0.35 Mean [{circumflex over (ν)}]_(x) [m/s] 0.067 0.2450.407 0.449 0.342 −0.303 −0.309 Std [{circumflex over (ν)}]_(x) [m/s]0.008 0.033 0.014 0.019 0.018 0.023 0.016 Error [%] — −2 −5.3 −10.2 2.213.4 11.7 Accuracy [%] 11.9 13.4 3.4 4.2 5.2 7.5 5.1 Values of the realtransverse velocity and the estimates mean and standard deviation for5000 scatterers and an emit focal distance of 150 mm. The error iscalculated as: (mean({circumflex over (ν)}_(x)) − v_(x))/v_(x), and theaccuracy is mean({circumflex over (ν)}_(x) − ν_(x))/std({circumflex over(ν)}_(x)).

TABLE 3 Parameters used for the parabolic profile simulation. PulseCenter Number of Sampling repetition Focal frequency Wavelength cyclesfrequency period distance f₀ [Hz] λ = c/f₀ [m] M f_(s) [Hz] T_(prf) [s]F [m] 3 10⁶ 0.51 · 10⁻³ 2 100 10⁶ 1/3500 150

TABLE 4 Values of the lateral and axial intervals used in thesimulations. Lateral line Lateral interval Axial length Axial intervallength [m] [m] [m] [m] 10 · 10⁻³ 0.01 · 10⁻³ Z(θ) 0.1 · 10⁻³

TABLE 5 Mean and standard deviation of the distance between the trueprofile and the estimated one. They are calculated for a differentcross- section distance for each of the angles so that only thesignificant part of the plots are taken into account: r ∈ [−8.5,6] for θ= 30°, r ∈ [−8.3,8.3] for θ = 45°, r ∈ [−8,8] for θ = 60°, r ∈ [−5,5]for θ = 90° (r indicates distance to the center of the vessel in [mm]).θ [deg] 30 45 60 90 Mean |{circumflex over (ν)} − ν| [m/s] 0.081 0.0560.04  0.073 Std |{circumflex over (ν)} − ν| [m/s] 0.047 0.038 0.0290.029

References

[1] D. J. Phillips, K. W. Beach, and J. Primozich D. E. Strandness.Should results of ultrasound Doppler studies be reported in units offrequency or velocity? Ultrasound Med. Biol., 15:205-212, 1989.

[2] U. Moser. Verfahren und Messanordnung zum Messen des Volumenstromesin einer Schicht mit reflektierender Struktur, European patentapplication, publication number 616 231 A2.

[3] M. D. Fox. Multiple crossed-beam ultrasound Doppler velocimetry.IEEE Trans. Son. Ultrason., SU-25:281-286, 1978.

[4] G. E. Trahey, J. W. Allison, and O. T. von Ramm. Angle independentultrasonic detection of blood flow. IEEE Trans. Biomed. Eng.,BME-34:965-967, 1987.

[5] V. L. Newhouse, D. Censor, T. Vontz, J. A. Cisneros, and B. B.Goldberg. Ultrasound Doppler probing of flows transverse with respect tobeam axis. IEEE Trans. Biomed. Eng., BME-34:779-788, 1987.

[6] O. Bonnefous. Measurement of the complete (3D) velocity vector ofblood flows. In Proc. IEEE Ultrason. Symp., pages 795-799, 1988.

[7] J. A. Jensen and N. B. Svendsen. Calculation of pressure fields fromarbitrarily shaped, apodized, and excited ultrasound transducers. IEEETrans. Ultrason., Ferroelec., Freq. Contr., 39:262-267, 1992.

[8] J. A. Jensen. Estimation of Blood Velocities Using Ultrasound: ASignal Processing Approach. Cambridge University Press, New York, 1996.

What is claimed is:
 1. An apparatus for measuring the velocity of amoving object or a collection of moving objects moving in apredetermined direction and at a predetermined distance, the apparatuscomprising: a generator for generating excitation signals, an emittingtransducer for transforming the excitation signals into wave energy andfor emitting said wave energy in a predetermined direction ofpropagation towards the object or objects, a receiving transducer forreceiving from said moving object or objects, signals generated byinteraction with said wave energy emitted from said emitting transducer,said emitting transducer and said receiving transducer having respectivesensitivities for selectively creating a plurality of focal points,which in combination define a focus line at the predetermined distanceand extending in the predetermined direction of movement, and means forcalculating a correction of signals received from the moving object orobjects on the focus line, and for calculating, based on thecorrelation, the velocity of the moving object or objects in thedirection of the focus line.
 2. An apparatus according to claim 1wherein said wave energy is pulsed wave.
 3. An apparatus according toclaim 1 wherein said wave energy is sound energy.
 4. An apparatusaccording to claim 1 wherein said sound energy is ultrasound energy. 5.An apparatus according to claim 1 wherein said wave energy iselectromagnetic energy.
 6. An apparatus according to claim 1 whereinsaid emitting transducer is an array transducer including a plurality ofemitting transducer elements.
 7. An apparatus according to claim 6further comprising an emit beam former for receiving said generatingexcitation signals and for supplying each of said plurality of emittingtransducer elements with individual excitation signals each having apredetermined time delay relative to the others of said individualexcitation signals.
 8. An apparatus according to claim 7 wherein saidindividual excitation signals have time delays resulting in focused waveenergy being emitted.
 9. An apparatus according to claim 1 wherein saidreceiving transducer is an array transducer including a plurality ofreceiving transducer elements.
 10. An apparatus according to claim 9further comprising a receive beam former for receiving signals from saidplurality of receiving transducer elements and for delaying each of saidsignals from said plurality of receiving transducer elementsindividually relative to the others of said signals from said pluralityof receiving transducer elements.
 11. A method for measuring thevelocity of a moving object or a collection of moving objects moving ina predetermined direction and at a predetermined distance, the methodcomprising: emitting excitation signals of pulses of wave energy in apredetermined direction of propagation towards the object or objects,receiving from said moving object or objects, reflected signalsresulting from interaction of emitted wave energy with the moving objector objects, the emission of excitation signals and the reception ofreflected signals having respective sensitivities for selectivelycreating a plurality of focal points, which in combination define afocus line at the predetermined distance and extending in thepredetermined direction, calculating a correlation of signals receivedfrom the moving object or objects on the focus line, and calculating,based on the correlation, the velocity of the moving object or objectsin the direction of the focus line.
 12. A method according to claim 11wherein said wave energy is pulsed wave energy.
 13. A method accordingto claim 11 wherein said wave energy is sound energy.
 14. A methodaccording to claim 11 wherein said sound energy is ultrasound energy.15. A method according to claim 11 wherein said wave energy iselectromagnetic energy.
 16. A method according to claim 11 wherein saidemitting transducer is an array transducer including a plurality ofemitting transducer elements.
 17. An method according to claim 16further comprising an emit beam former for receiving said generatingexcitation signals and for supplying each of said plurality of emittingtransducer elements with individual excitation signals each having apredetermined time delay relative to the others of said individualexcitation signals.
 18. A method according to claim 17 wherein saidindividual excitation signals have time delays resulting in focused waveenergy being emitted.
 19. A method according to claim 11 wherein saidreceiving transducer is an array transducer including a plurality ofreceiving transducer elements.
 20. A method according to claim 19further comprising a receive beam former for receiving signals from saidplurality of receiving transducer elements and for delaying each of saidsignals from said plurality of receiving transducer elementsindividually relative to the others of said signals from said pluralityof receiving transducer elements.